/*
 * complex_li.c
 * Copyright (C) 2014 Kevin Marshall Stueve
 *
 * This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 * Author: Kevin Marshall Stueve; August 2014
 */
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "complex_li.h"

/*
 * Code to import the zeta zeros to an array
 *     Author: Kevin Marshall Stueve; August 2014
 *     Source: http://www.cplusplus.com/reference/cstdio/fgets/
 *
 */
int import_zeros() {
	const char *filename ="c:\\prime_data\\zeros6";
	FILE *file = fopen(filename, "r");
	if (file ==NULL) perror("Error opening file");
	char line [100];
	int line_num=0;

	while(!feof(file)){
        if (fgets(line,100,file)!=NULL) {
        	zeta_zeros[line_num] = atof(line);
        	line_num += 1;
        }
	}
	fclose(file);
	return line_num;
}

/*
 * pari-gp computation of li(x^{1/2+I t})+li(x^{1/2-I t})
 *     Author: Tomas Oliveira e Silva January 2010
 *
 * li(x^\rho) = x^\rho/u*(1+1/u+2!/u^2+3!/u^3+...),
 * with \rho=1/2+it and u=\rho\log(x)
 *
 *	li(x,t)=
 *	{ local(u,s0,s1,s2,lx,tol,k);
 *	  u=(1/2+I*t)*log(x); //log(x^{1/2+I t})
 *	  s0=2*x^(1/2+I*t)/u;
 *	  tol=1e-9/abs(s0); //we desire an error of less than 10^{-9}
 *	  s1=s2=1;
 *	  k=1;
 *	  while(1,
 *	    if(k>100,
 *	      error("li: unable to attain the desired precision");
 *	    );
 *	    s2*=k/u;
 *      s1+=s2;
 *      if(abs(s2)<tol,
 *	      break();
 *	    );
 *	    k++;
 *    );
 *	  return(real(s0*s1));
 *	}
 */

/*
 * C optimized calculation of li(x^(0.5+it))+li(x^(0.5-it))
 *     Author: Kevin Stueve; March 2010
 */
long double li_rho(long double x, long double t) {
	long double logx = log(x);
	long double sqrtx = sqrt(x);
	long double u_real = 0.5 * logx;
	long double u_imag = t * logx;
	long double u_norm = u_real * u_real + u_imag * u_imag;
	long double so_real_num = 2 * sqrtx * cos(t * logx);
	long double so_imag_num = 2 * sqrtx * sin(t * logx);
	long double so_real = (so_real_num * u_real + so_imag_num * u_imag) / u_norm;
	long double so_imag = (so_imag_num * u_real - so_real_num * u_imag) / u_norm;
	long double tol=pow(10,-9)/(sqrt(so_real * so_real + so_imag * so_imag));
	long double tol_squared=tol*tol;
	long double s1_real = 1;
	long double s1_imag = 0;
	long double s2_real = 1;
	long double s2_imag = 0;
	long double s2_real_temp;
	long double s2_imag_temp;
	int k = 1;
	long double norm;
	for(;;) {
		if (k>100) {
			printf("Error on t=%Lf\n",t);
			return 0;
		}
		s2_real_temp = k * (s2_real * u_real + s2_imag * u_imag) / u_norm;
		s2_imag_temp = k * (s2_imag * u_real - s2_real * u_imag) / u_norm;
		s2_real = s2_real_temp;
		s2_imag = s2_imag_temp;
		s1_real += s2_real;
		s1_imag += s2_imag;
		norm = s2_real * s2_real + s2_imag * s2_imag;
		if (norm < tol_squared) {
			break;
		}
		k++;
	}
	return so_real * s1_real - so_imag * s1_imag;
}

 /* C optimized calculation of li(x^(1.5+it))+li(x^(1.5-it))
  * li(x^(1.5+it))=li(x^(rho+1))
  * See code and comments for li_rho above.
  *     Author: Kevin Marshall Stueve; August 2014
  */
 long double li_rho_plus_one(long double x, long double t) {
 	long double logx = log(x);
 	long double sqrtx = sqrt(x);
 	long double u_real = 1.5 * logx;
 	long double u_imag = t * logx;
 	long double u_norm = u_real * u_real + u_imag * u_imag;
 	long double so_real_num = 2 * x * sqrtx * cos(t * logx);
 	long double so_imag_num = 2 * x * sqrtx * sin(t * logx);
 	long double so_real = (so_real_num * u_real + so_imag_num * u_imag) / u_norm;
 	long double so_imag = (so_imag_num * u_real - so_real_num * u_imag) / u_norm;
 	long double tol=pow(10,-9)/(sqrt(so_real * so_real + so_imag * so_imag));
 	long double tol_squared=tol*tol;
 	long double s1_real = 1;
 	long double s1_imag = 0;
 	long double s2_real = 1;
 	long double s2_imag = 0;
 	long double s2_real_temp;
 	long double s2_imag_temp;
 	int k = 1;
 	long double norm;
 	for(;;) {
 		if (k>100) {
 			printf("Error on t=%Lf\n",t);
 			return 0;
 		}
 		s2_real_temp = k * (s2_real * u_real + s2_imag * u_imag) / u_norm;
 		s2_imag_temp = k * (s2_imag * u_real - s2_real * u_imag) / u_norm;
 		s2_real = s2_real_temp;
 		s2_imag = s2_imag_temp;
 		s1_real += s2_real;
 		s1_imag += s2_imag;
 		norm = s2_real * s2_real + s2_imag * s2_imag;
 		if (norm < tol_squared) {
 			break;
 		}
 		k++;
 	}
 	return so_real * s1_real - so_imag * s1_imag;
 }

 /* C optimized calculation of x*li(x^(0.5+it))+x*li(x^(0.5-it))
  * See code and comments for li_rho above.
  *     Author: Kevin Marshall Stueve; August 2014
  */
 long double x_li_rho(long double x, long double t) {
  	long double logx = log(x);
  	long double sqrtx = sqrt(x);
  	long double u_real = 0.5 * logx;
  	long double u_imag = t * logx;
  	long double u_norm = u_real * u_real + u_imag * u_imag;
  	long double so_real_num = 2 * x * sqrtx * cos(t * logx);
  	long double so_imag_num = 2 * x * sqrtx * sin(t * logx);
  	long double so_real = (so_real_num * u_real + so_imag_num * u_imag) / u_norm;
  	long double so_imag = (so_imag_num * u_real - so_real_num * u_imag) / u_norm;
  	long double tol=pow(10,-9)/(sqrt(so_real * so_real + so_imag * so_imag));
  	long double tol_squared=tol*tol;
  	long double s1_real = 1;
  	long double s1_imag = 0;
  	long double s2_real = 1;
  	long double s2_imag = 0;
  	long double s2_real_temp;
  	long double s2_imag_temp;
  	int k = 1;
  	double long norm;
  	for(;;) {
  		if (k>100) {
  			printf("Error on t=%Lf\n",t);
  			return 0;
  		}
  		s2_real_temp = k * (s2_real * u_real + s2_imag * u_imag) / u_norm;
  		s2_imag_temp = k * (s2_imag * u_real - s2_real * u_imag) / u_norm;
  		s2_real = s2_real_temp; s2_imag = s2_imag_temp;
  		s1_real += s2_real;
  		s1_imag += s2_imag;
  		norm=s2_real * s2_real + s2_imag * s2_imag;
  		if (norm < tol_squared) {
  			break;
  		}
  		k++;
  	}
	return so_real * s1_real - so_imag * s1_imag;
  }
